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ity and their ability to reproduce various spiking patterns of cortical neurons, and 



O 



X 



are particularly used for large network simulations. These models describe the dy- 
namics of the membrane potential by a nonlinear differential equation that blows 
up in finite time, coupled to a second equation for adaptation. Spikes are emitted 
when the membrane potential blows up or reaches a cutoff 6. The precise simula- 
tion of the spike times and of the adaptation variable is critical for it governs the 
spike pattern produced, and is hard to compute accurately because of the explod- 
ing nature of the system at the spike times. We thoroughly study the precision of 
fixed time-step integration schemes for this type of models and demonstrate that 
these methods produce systematic errors that are unbounded, as the cutoff value is 
increased, in the evaluation of the two crucial quantities: the spike time and the 
value of the adaptation variable at this time. Precise evaluation of these quanti- 
ties therefore involve very small time steps and long simulation times. In order to 



achieve a fixed absolute precision in a reasonable computational time, we propose 
here a new algorithm to simulate these systems based on a variable integration 
step method that either integrates the original ordinary differential equation or the 
equation of the orbits in the phase plane, and compare this algorithm with fixed 
time-step Euler scheme and other more accurate simulation algorithms. 



1 Introduction 



The increasingly finer description of the biophysics of neurons and their ionic chan- 
nels now makes the biological mechanism of action potential firing understood in 
great detail^] This better understanding allows to propose very precise models 
of neurons mimicking the dynamics of ionic transfers through the channels. Yet, 



simple neuron models such as the integrate-and-fire model (Lapicque, 1907; Ger 



stner and Kistler, 2002) remain very popular in the computational neuroscience 



community, because they can be simulated very efficiently and, perhaps more im- 
portantly, because they are easier to understand and analyze. The drawback is 
that these simple models cannot account for the variety of electrophysiological 



behaviors of real neurons (see e.g. (Markram et al. , 2004) for interneurons) . Re 



cently, several authors introduced two- variable spiking models (Izhikevich, 2004 



Brette and Gerstner, 2005 Touboul, 2008a) which, despite their simplicity, can 



reproduce a variety of electrophysiological signatures such as bursting or regular 
spiking. Different sets of parameter values correspond to different electrophysio- 
logical classes. 

These models seem to provide a compromise between simplicity and versa- 



tility. Simplicity, in the sense that it allows analytical studies (Touboul, 2008a 



Touboul and Brette, 2009), efficient numerical simulations for very large networks 



simulations (Izhikevich and Edelman, 2008), and a small number of parameters. 



And versatility, because it is able to reproduce a variety of possible behaviors that 



neurons present: in (Izhikevich, 2004 Touboul, 2008a), the authors provide simu- 



lations corresponding to a wide variety of neuronal behaviors, and simple sets of 
parameters corresponding to each behaviors can be determined, either analytically 
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Many questions still remain open, as for instance reviewed in (Hille 



2001). 
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flTouboull |2008a| iTouboul and Brettel 120091) or numerically flNaud et al.l 120081) 



and that can be precisely and easily tuned to fit intracellular recordings ( Clopath 



et al. 


2007 


Jolivet et al. 


2008alb 


Cyrille et al. 


2010) 



the interest of mathematicians, for the nonlinear nature of the dynamics combined 
with the discrete nature of spikes makes of this class of models an interesting novel 



mathematical object (Touboul, 2008a; Touboul and Brette, 2009) 



These models describe the dynamics of the membrane potential by a nonlinear 
differential equation, coupled to a second equation for adaptation, and a separate 
discrete mechanism accounts for the spike emission. In details, these models satisfy 
reduced equations of the type: 



dv 
dt 

dw 
dt 



F(v) - w + I 
a(bv — w) 



(1) 



where a and b are non-negative parameters accounting respectively for the time 
scale of the adaptation variable with respect to the membrane potential's time 
scale and for the coupling strength between the two variables. The function F(-) 
is a real function, that satisfies the following assumptions: 

Assumption (Al). 

1. F is regular (at least three times continuously differentiable) 

2. F is strictly convex, and lim F\x) < 0, 

x— >— oo 

3. There exists e > and a > for which F(v) > av 1+£ when v — > oo (we will 
say that F grows faster than v 1+e when v — y oo). 

Under these assumptions, the membrane potential blows up in finite time for 
some initial conditions (see (Touboul and Brette 2009) and appendix [A]) . A spike 
is emitted at the time t* when the membrane potential v reaches a cutoff value 9 
or when it blows up. At this time, the membrane potential is reset to a constant 
value c and the adaptation variable is updated to w(t*)+d where w(t*) is the value 
of the adaptation variable at the time of the spike and d > is the spike-triggered 
adaptation parameter. Furthermore, if the function F satisfies the assumption: 
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Assumption (A2). There exists e > and a > for which F(v) > av 2+e when 
v — > oo (F grows faster than v 2+£ when v — > oo), 

then the adaptation variables converges to a finite value when the membrane 
potential blows up, which allows replacement of the strict voltage threshold of 



classical (linear) integrate-and-fire neurons (Lapicque, 1907 Stein, 1967) by a more 



realistic smooth spike initiation (see e.g. (Latham et al. 2000 Fourcaud-Trocme 



et al. 2003)). 



Among these models, the widely used quadratic adaptive model (Izhikevich 
2004) corresponds to the case where F(v) = v 2 , and has been recently used by 



Eugene Izhikevich and coworkers (Izhikevich and Edelman, 2008) in very large 



scale simulations of neural networks. This model does not satisfies assumption 
(A2)| and in that case the adaptation variable w blows up at the times of the spikes, 
which implies that the spike patterns produced by the simulation of this model 
sensitively depends on the choice of the cutoff, parameter which does not have any 
biological interpretation (see (Touboul, 2009[ )). The adaptive exponential model 
(Brette and Gerstner, 2005) corresponds to the case where F(v) = e v — v. It has 
the interest that its parameters can be related to electrophysiological quantities, 



and has been successfully fit to intracellular recordings of pyramidal cells ( Clopath 



2007 


Jolivet et al. 


2008a 


Badel et al., 2008) 



2008b) corresponds to the case where F(v) — v + olv . It has the advantage 
of being able to reproduce all the behaviors featured by the other two and also 
self-sustained subthreshold oscillations which are of particular interest to model 
certain nerve cells. 

In these models, the neural code is assumed to be contained in the times 
of the spikes. Therefore, computing the spike times with accuracy is essential. 
Moreover, the reset mechanism makes critical the value of the adaptation variable 
at the time of the spike. Indeed, when a spike is emitted at time t*, the new initial 
condition of the system ([!]) is (c, w(t*) + d). Therefore, this value w(t*) governs 
the subsequent evolution of the membrane potential, and hence the spike pattern 



produced. For instance in (Touboul and Brette, 2008 Touboul, 2009; Touboul 



and Brette, 2009), the authors show that the sequence of reset locations after each 



spike time shapes the spiking signature of the neuron. They also show that small 
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perturbations can result in dramatic changes, since the spike pattern produced 
undergoes bifurcations as a function of the parameters. This discontinuity linked 
with the spike emission might lead to error accumulation on a spike train emission. 

The explosion of the membrane potential variable (and possibly of the adap- 
tation variable) makes the simulation of these models very delicate. Indeed, the 
explosion implies that the membrane potential v , as a function of time, presents 
singularities at the spike times, where in particular both the membrane potential 
and the adaptation variables have an infinite derivative (the explosion in finite 

time always corresponds to an infinite derivative, and equation ([I]) implies that 

[9l 

at these times the derivative of w also tends to infinity ). This property of the 
model makes the accurate simulation of these models very difficult. For one dimen- 
sional models where no adaptation (variable w) is taken into account, the spike 



time can be corrected following (Fourcaud-Trocme et al. 2003) by computing the 
explosion time in closed form, which allows a good evaluation of the spike times. 
Unfortunately, such corrections do not exist for the more biologically realistic bidi- 
mensional models with adaptation, and the precise evaluation of spike times and 
adaptation value at these times necessitate the development a new computational 
techniques. 

In this paper, we start by thoroughly studying the precision of fixed time 
step methods such as the Euler integration scheme that are the most widely used 
methods to simulate such equations. We will observe that the error made is 
unbounded as the threshold 6 is increased, which leads us to propose in section 
[3] a new algorithm for accurately and efficiently simulating the solutions of these 
equations, and in particular two crucial features of the model, namely the spike 
time and the value of the adaptation variable at these times. The comparison of 
this new algorithm with some other algorithms is then discussed. All along the 
paper, we will make use of some theoretical results on the orbits of the dynamical 
system ([I]) that are summarized and proved in appendix |Aj 



o 

^it can be easily proved that w(t) is always negligible compared to v(t) (i.e. w(t) = o(v(t) j) 
at the explosion times 
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2 Fixed Time Step Simulation Methods 



Most simulation algorithms found in the literature for simulating neuron mod- 
els of the type of equation Q involve direct integration of the equations with a 



fixed time step integration algorithms such as the Euler scheme (Izhikevich, 2003 



Clopath et al., 2008 Touboul and Brette, 2008), or more precise Runge-Kutta 



like schemes. In these simulation algorithms, spikes are emitted when the voltage 
variable reaches a given threshold 9. The use of such fixed time-step methods for 
solving this type of blowing-up equations was initially motivated by the observa- 
tion that the spike time is easily evaluated because of the explosion property of the 



membrane potential variable v, as discussed in (Izhikevich, 2007), where a typical 
Euler algorithm is suggested. However, recent research established that in order 
to accurately simulate these models, one needs to also have an important precision 
in the evaluation of the adaptation variable at the time of the spike. Indeed, it 



was shown that this variable governs spiking patterns (Touboul, 2009 Touboul 



and Brette, 2009), and that small changes in the evaluation of this variable can 



quantitatively and qualitatively modify them. 

The crucial question that this observation raises is that beyond the necessary 
precision in the spike time, the model requires to accurately evaluate the value of 
the simulated adaptation variable at the time of the spike. In this view, though 
usual fixed time-step methods might provide a fair precision on the evaluation 
of the spike time (we will show in the sequel that even this evaluation is very 
imprecise with such fixed time-step methods, involving unbounded errors as the 
spiking cutoff is increased), they are very imprecise for calculating the value of 
the adaptation variable at the time of the spike as we rigorously show in the 
present section. This imprecision will substantially modify the value of the time 
of the next spikes, and this effect will be increasingly important as the number of 
spikes emitted increases. Indeed, because of the discontinuity induced by the spike 
emission and the reset, the numerical errors will accumulate. We first analyse the 
precision for the computation of the first spike, before addressing the question of 
how do these errors accumulate. 
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2.1 Precision of the computation of the first spike time 
and adaptation variable at this time 

We now go into the mathematical details of the precision analysis of such fixed 
time step methods for simulating the solutions of equation 0. We concentrate 
on the simpler and most widely used fixed-time step Euler scheme, with threshold 
9 and time step r. The numerical solutions computed through such numerical 
procedure are solutions of the recursion: 

v n +i = v n + r(F(v n ) -w n + I{nr)) 

(2) 

w n+1 = w n + ra(bv n -w n ) 

starting form a given initial condition (v ,w ). We denote X n := (v n ,w n ) and $ 
the map such that X n+ \ = given by the recursion @. We define a zone of 



the phase plane Z* as a subset of the spiking zone defined in (Touboul and Brette 



2009) which will be of particular interest in the sequel. 



Definition 1. Let Z(V) be the region of the phase plane (v,w) defined as: 

Ziy) = {(v, w) G R 2 ; v > V and w < bv}. 

Note that Z(V) constitute a decreasing family of nested sets. 

We further define, when the input current I(t) > I* is lowerbounded, 

Z*:= Z(r,b)=U v > v+{I * ib) Z(V), 

where v + (I*, b) is the largest fixed point of the system ([!]) when it exists, and — oo 
if it does not (see appendix [A]). We have if m(b) denotes the minimal value of 
F(v)-bv: 

• If I* > -m(fo), then Z(I*,b) = {(v,w) G R 2 ; w < bv}. 

• If I* < -m(b), then Z(I*,b) = {(v,w) G R 2 ; v > v + {I*,b) and w < bv} 

When one does not consider any threshold on the dynamics of the membrane 
potential variable v, any spiking trajectory will eventually belong to Z* after a 



transient phase. 

Lemma 1. Let us assume that (vq,wo) G Z* and ra < 1. Then for all n > 0, we 

have X n G Z* and both sequences (v n ) and (w n ) are strictly increasing and finite. 
Moreover, Hindoo v n = oo. 

Proof. We reason by induction on n. We know that X = (v ,w ) G Z*. Let 
us assume that for some n > 0, we have X n = (v n ,w n ) G Z*. The finiteness of 
X n+ \ = $(X n ) stems straightforwardly from the finiteness of X n . 

Moreover, since X n G Z* implies that w n < bv n and v n > v + (I*,b) when it 
exists, we have: 

F(v n ) -w n + I(nr) > F(v n ) - bv n + I* > 

and because of the convexity assumption made of F we have F(v) — bv + I* > 
for all v > v + (I*,b) when it exists, or for all v when no fixed point exists. This 
directly implies that v n+ i > v n . Moreover, since in Z*, we have w < bv, it is 
immediate to show that w n+ i > w n , implying that both sequences are increasing. 

We therefore have v n+ i > v n > v > v + (I*, b). Therefore, showing that X n+1 = 
$(X n ) G Z* only amounts proving that w n+ \ < bv n+ \. We have: 

w n+ i - bv n+1 = (w n - bv n ) + r[a (bv n - w n ) -b(F(v n ) - w n + I(nr))] 




>o <o 



< 0. 

Therefore, for all n > we have X n G Z*. 

To end the proof of the lemma, we therefore need to prove the unboundedness 
property of v n . Since we have i>„+i — v n > r{F{v n ) — bv n + /*) > r(F(v ) — bvo + 
/*) =: K wich is strictly positive constant, we have for all p G N the inequality 
v p > v + pK , and therefore for all M > there exists n G N such that v n > M 
which ends the proof of the lemma. □ 

We therefore observe that for any initial condition in the spiking zone, the 
adaptation variable of the solution of the Euler simulation scheme diverges (i.e. 
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v n — y oo when n — > oo). However, for any n G N, we have shown that t>„ is 
finite, hence the Euler scheme does not blow up in finite time. This fact implies 
in particular that the error made on the evaluation of v is a priori unbounded 
if there is no cutoff in the simulation: solutions of ([!]) starting from some initial 
condition (vq, wq) G Z* blow up in finite time (say t*) and the approximated value 
of v at time t* computed from the Euler scheme is finite, producing an infinite 
error on the evaluation of v. 

However, as we already mentioned, numerical methods truncate the solution 
as soon at the value of the adaptation variable reaches a given threshold (or 
cutoff) 9. The spike time corresponds to the first time the simulated trajectory 
exceeds 9, and is then reseted. For v < 9, the vector field Lipschitz-continuous, 
and it is known from standard numerical analysis theory that the thresholded 
discretization scheme is of order one in r, meaning that the error is bounded by 
a constant depending on the parameters of the model multiplied by r. Ensuring 
an absolute given precision on the spike time and the adaptation variable at this 
time therefore amounts controlling the value of the constant involved in the error, 
which is what we now do. 

Let us now have a closer look at the absolute error made in the thresholded 
discretized solution. We define the functions a(-) and /?(•) corresponding to the 
first order error in r of the discretization of v and w as: 

v n = v(nr) + t a(n r) + 0(r 2 ) 
w n = w(nr) + Tf3(n r) + 0(r 2 ) 

We have: 

a(0) = /3(0) = 0, 



and moreover (see e.g. (Hairer and Lubich, 1984 Isaacson and Keller, 1994 Sanz 



Serna and Verwer 



1986)): 



v n +i — v n = v(nr + r) — v(nr) + r(a(nr + r) — a{nr)) + 0(r 2 ) 
= rv'inr) H v"(tit) + T 2 a'(nr) + 0(r 3 ) 



2 



r 2 



r (F(w(nr)) - u>(nr) + /(nr)) + -/(nr) + rV(nr) + 0(r 3 ). 

By the recursion relationship between w„ and t> n +i we have: 

v n +i -v n = r {F(v n ) -w n + I{nr)) 

= t (F(v(nr) + ra(nr) + 0(r 2 )) - w{nr) - t(5(ut) + I(nr) + 0(r 2 )) 
= r (nr)) + F'(v(nr))ra(nT) - w(nr) - r/3(nr) + I(nr)) + 0(r 3 ) 

Equalizing the two expressions, we get: 

a'{nr) = F'(v(nr)) a(nr) — (3(nr) — -v"(tit) 

proceeding in the exact same fashion for the difference w n+ \ — w n , we are led to 
the following system of ordinary differential equations on the first order error of 
the Euler scheme: 

' a'it) =F\v(t))a(t)-P(t)-\v"(t) 
f3'(t) = a (bait)- (3(t))-\w"(t) 

Using the equations ([TJ governing v and w, we obtain the equations: 



(4) 



a' = F'(v) a - (3 - \\F\v) (F(v) - w + I) - a (bv - w) + /' 
P' = a (ba-p)- %\b(F(v)-w + I)-a(bv-w) 



It is easy to prove using Gronwall's theorem along the same lines as in (Touboul 



2009 ) that both errors tend to infinity when no threshold is considered, as we also 
noted from the fact that the continuous solution v (t) blows up in finite time and 
the discretized version v n remains finite for all n G N. 

Let us consider a spiking trajectory. From theorem [2] of appendix [A] and in 
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lemma [TJ we know that as soon as the trajectory enters the spiking zone Z* , 
the error can be parametrized as a function of the membrane potential v(t). In 
the sequel, instead of studying the error functions a(t) and (3(t), it will appear 
particularly convenient to study the composed applications A(v) = (a o T)(v) = 
a(T(v)) and B(v) = (/3 o T)(y) = (3{T{v)) where T(v) is defined in theorem [2] of 



appendix A as the inverse of the function 1 1- >■ v(t) and o denotes the composition 



of applications. These two functions satisfy the equations: 



dA ___ p 

dv F( 



(F(v)-W + IoT)-a(bv-W)+I'oT] 
"- } (bA-B)-^{b (F(v) - W + IoT)-a(bv-W)} 



dB _ 
dv F( v ) 

These equations, similarly to initial model's equations, are quite hard to solve in 
their general setting. However, using the fact that W < bv and that on a spiking 
trajectory W = o(v) and similarly j3 = o(v), the behavior of the error close to a 
relatively high threshold can be approximated by the solution of the equations: 



dA _ F'(v) a _ 1 pi 
dv F(v) 2 r 



(5) 

(h A — U\ — In h 

dv F(v) 



- ^{bA-B)-\ab 



These equations can be easily understood heuristically. Indeed, since we known 
that at the times of the spikes, the membrane potential blows up and the adap- 
tation variable is dominated by the value of the membrane potential variable, the 
error made on the membrane potential evaluation essentially stems from the di- 
vergence membrane potential variable, that asymptotically satisfies an equation 
of type y = F(y) + /. This is why the first equation is exactly the same as the one 
for the error of fixed-time step methods for this one-dimensional equation. The 
error A computed is the term that mostly disturbs the evaluation of the adap- 
tation variable, and acts on the second equation of (|5| as an independent input 
function. This is coherent with the heuristic argument that the main part of the 
errors made on the estimation of the adaptation variable at the time of the spike 
is therefore related to the imprecision in the evaluation of the membrane potential 
variable close to the explosion. 

The study of the first equation of ^ and of the threshold crossing time is 
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quite straightforward at this point. The first equation of (J5J) is integrated as: 
A(v) = i(log (F(v )) - \og(F(v))) F(v) 



which yields for the error on the adaptation variable: 



B(v) 




e i: nk) dui (l og (F(u)) - \og{F{v )))du 



and the error made on the evaluation of the adaptation variable because of the 
discretization at the threshold crossing time is therefore equal to \B{9)\. It is now 
easy to instantiate the models and find an approximation of the error: 

(i). For F(v) having a polynomial dominant term v m (which is for instance the 
case of Izhikevich quadratic model and the quartic model), we have: 



(ii). For F(v) equivalent to an exponential function (as it is the case in the 
adaptive exponential model), we have the error estimate: 



bounded and diverging as 6 increases. This result implies that the greater the 
cutoff is the more imprecise the evaluation of the adaptation variable at the cutoff 
time is, and that the error diverges as the cutoff value is increased. Hence the 
use of fixed time step methods to approximate such nonlinear integrate-and-fire 
models leads to substantial quantitative errors on the evaluation of v and w that 
might change the qualitative behaviors, as we will further address in section |4j 
We also note that the error made is always negative: the evaluated w at the time 
of the spike will be systematically smaller than the actual value of the adaptation 




) 



du 




These functions are plotted as a function of v in Figure [T] They are both un 
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Figure 1: Error made on the estimation of the adaptation variable as a function of 
the threshold used. We see a fast divergence of the error. Red solid curve: quartic 
case, Green dotted curve: exponential case 



We now adress the following question: is the spike time accurately evaluated 
by this procedure? The answer to this question was understood to be positive, 
and the argument given was the divergence of the membrane potential at the 
time of the spike. However, as we rigorously prove in the sequel, this is not 
exact. Indeed, the explosion only concerns the differential equation, whereas the 
numerical procedure does not blows up in finite time, but diverges when time 
tends to infinity. This induces increasingly large (and even infinite) errors on the 
evaluated membrane potential v that produces systematic errors at the threshold 
crossing time, and these errors are unbounded as the cutoff is increased. The proof 
of this fact is quite difficult to produce in the full system. However, as already 
discussed, because of the fact that during the explosion we have w = o(v), we will 
restrict the study to the comparison of the recursion y n +i = y n + F{y n ) and the 
blowing up solution of y' = F(y), and only for power and exponential F functions, 
which cover most practical cases. 

Let us treat in parallel the power and the exponential cases. Let m > 1 
be a real number, and consider respectively the ordinary differential equations: 



variable. 
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The solution of these ordinary differential equation reads (respectively): 



y P (t)=Hr m -(m-l)t] 1/(1 - m) 
y e (t)=y e -\og(l-te^) 

and blow up at time t* = {yo) l m /{ m — 1) > (resp t* = e~ y o). The numerical 
solution from Euler discretization satisfy the recurrence equation: 

L p n+1 =yl + r{yir 

[Vn+i =y e n + rey- 

with initial condition y% and y% respectively The approximation error produced 
by the discretization algorithm was previously computed, and are equal to: 

AP{t) = _a 1/ P( t )m log ^^ 
A e(t) = -±(s/ c (*)-J/§)e* e W 

for < t < t* , and therefore we have: 

yl = yP(t) - rf yP(t) m log (^) + 0(r 2 ) 
yl = y e (t) - r\ (y e (t) - yl) e ™ + 0(r 2 ) 

when r — > and for < t < t*. Let us now define z p (k) (resp. z e (k)) the 
continuous variable linear interpolation of y? (resp y^) between two integers n 
(we have y e {n) = z e {n) for all n G N and between two consecutive integers the 
function k G [n, n + 1] h- >■ z e (k) is linear, and similarly for the polynomial case 
with indexes p). The relation between z p (kr) and y p (t) (resp. z e (kr) and y e (t)) 
is the same as the relation between y p (t) and y p (resp. y e (t) and y^) since the 
errors made by a linear interpolation in an interval of length r are of order 0{t 2 ). 
Let us now define n* G R and n* G R such that z e (n* e ) = z p (n* p ) = 9. We aim at 
comparing n* and n** defined by y e {n**r) = 9 (and similarly in the polynomial 
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case). We have: 

= (yg) 1 -" _ i 

h V (m-l)r (m-l)0 m - 1 T 

n* e * = i (e"»S - e - e ) 
Let now Z u = y u (n*r) (u = e,p), we have: 

= Z p -|rmZ™log(§) 
=Z e -l(Z e -y e )e z * 

which can be solved up to the second order (by identifying the constant and the 
linear in r coefficients of the solution): 

z p =e + ^fe m \o g (^) + o(r 2 ) 

Ze =e + r(e-y^e e + 0(r 2 ) 
and eventually inverting the expression of y(t) one obtains: 

K =n;* + f^log(^) 

Therefore, there is a systematic delay of the estimated spike time produced com- 
pared to the actual spike time, that increases as the cutoff is increased. In order 
to achieve precision on the evaluation of this variable also, one therefore needs to 
use small time steps. 

Therefore, we showed, by analogy to one dimensional ordinary differential equa- 
tion that capture the explosion of the membrane potential at the time of the spike, 
that a systematic error is made on the evaluation of the spike time, which is of first 
order in r with a bounded coefficient that diverges when the cutoff is increased, 
and most important that the errors made in the estimation of the adaptation 
variable at the time of the spike, of order one in r also, have a coefficient that is 
of large amplitude and that is fast diverging when the cutoff is increased. This 
implies that very small time steps are necessary to achieve a given precision of the 
simulation, which leads to excessive precision in region of parameters where the 
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vector field is smooth and that questions the efficiency of the obtained simulation 
algorithm. 

We observe that the faster the divergence of the membrane potential, linked 
with the growth of the function F(v), the larger the error made on the adaptation 
variable at the cutoff time, and the smaller the error made on the spike time are. 

2.2 Error Accumulation during a spike train 

We have seen that starting from the same initial condition, a fixed time step Euler 
scheme produces errors on the evaluation of the spike time and on the value of 
the adaptation variable at this time. Therefore, when simulating the emission of 
a spike train, an additional deviation between the numerical spike train computed 
and the actual spike train appears linked with a shift in the spike time and on the 
value of the adaptation variable at these times, and as we show here, since the 
map is expanding, the distance between two trajectories increases as time goes 
by, which increases numerical errors produced. Indeed, assume that the actual 
system spikes at time t* and the simulated system at time £**, and that the actual 
value of the adaptation at this time w* is approximated by w**. We therefore 
need to compare the solutions of equation (JT1) with different initial conditions 
(t*,c,w*), solution we denote (v 1 ,wi), and (t**, c, to**), denoted (^2,^2)- We have 
for t > max(t*,t**): 



The error on the value of the adaptation variable w is therefore governed by the 
error made on the membrane potential variable, and this difference is unboundedly 
large close to the spike emission time where the flow is very expanding because 
of assumption |(Al)| iii). Indeed, for all K > there exists B > such that for 
all x,y greater than B, \f{x) — f(y)\ > K\x — y\). Therefore the errors on the 
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computation of the time and adaptation variable are expanded if the threshold is 
large enough, and there is no way to control it since this divergence is linked with 
the divergence of the solutions of equation Q with different initial conditions 

We therefore showed in this section that the Euler fixed time step methods 
yielded important errors on the computation of the spike times and on the value 
of the adaptation variable at these times. The results still hold for any fixed time 
step algorithm since these are linked with the very nature of the problem, namely 
the explosion at the time of the spike and the reset discontinuity, in addition to the 
infinite slope of the adaptation variable at the spike time. These issues motivate 
the introduction of an integration scheme adapted to the nature of the equations 
we simulate. 



3 Fixed Absolute Precision Simulation Algorithm 

The explosion, and the infinite slope of the adaptation variable at the time of the 
spike make the accurate simulation of the spike very delicate, and usual fixed time- 
step methods necessitate small time steps in order to achieve a precise evaluation 
of the spike time and of the adaptation variable at these times, which results in 
excessive precision in regions where the vector field varies slowly. 

In this section we propose a novel algorithm in order to overcome these difficul- 
ties. The heuristic idea behind this approach is the following: in order to overcome 
the difficulty of computing close to the explosion the diverging membrane poten- 
tial variable as a function of time, we propose to turn the picture around and to 
simulate time as a function of the membrane potential. The simulation algorithm 
we propose is based on a simulation of the trajectories in the phase plane as stud- 
ied in appendix |Aj in particular in theorem [2] More specifically, the algorithm is 
based on the local inversion theorem. In details, let M > be a fixed constant, 
to G R a given time and assume that (v(t),w(t)) is an orbit of ([TJ containing 
the point (to,v ,w ). Provided that v'(t ) = F(v ) — w + /(to) > M, there ex- 
ists 5 > such that t i— >■ v(t) is invertible on [t — 5,t + 5}. Let us denote by 



°Note also that this effect is increasingly important as the cutoff is increased. The complete 
proof of this result involves Gronwall theorem, and can be related to the contraction theory (see 
e-g 



Slotine and Li (19911; Lohmiller and Slotine (1998)) 
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v i — y T(v) the inverse of the function t i— > v(t) and by W(v) the composed appli- 
cation w o T(v) = w(T(v)). These functions are differentiable in a neighborhood 
of (vo,Wq) and the differential reads: 

'dT 1 
cb ~ - W + I(T) 

(6) 

a(bv-W) 
,^kT ~ - W + I{T) 

This inversion is valid as long as v'(t) is invertible, i.e. as long as F(v(t)) —w(t) + 
I{t) > 0. We have 

v "(t) = F'{v{t)){F{v{t)) - w{t) + I(t)) - a {bv{t) - w(t)) + l\t) 

and therefore 

\v"(t)\ < \F'(v)(F(v) - W(v) + I(T(v))) -a(bv- W(v)) + I'(T(v)) 

with F(vq) — w + /(to) > We can therefore define 5' > as a function of this 
second derivative, and such that the description (|6| is valid on [vq — 8', Vq + 5'}. 
This choice of 5' will be further discussed in the description of the algorithm. 

Moreover, it is important to note that as soon as the trajectory enters the 
spiking zone, in particular when it enters the zone Z*, this description of the 
trajectories in the phase plane is valid until the spike emission (explosion of the 
membrane potential). More precisely, we demonstrate in appendix|A]the following: 

Theorem 2. Assume that the input current I(t) depends on time, and moreover 
that t i — y I(t) is lower bounded (at least on the time interval considered {tq,t\)), 
i.e. I(t) > I* for any t G (to,Ti). Then we have: 

(i) . Z(I*,b) =: Z* is stable under the flow of the equation 

(ii) . Let (vo,wq) E Z* , and denote (v(t),w(t)) the solution of equations ([TJ having 

initial condition (v ,w ) at time t G (r ,ri). Then t \-t v{t) is strictly 
increasing and has a strictly positive derivative, and therefore is smoothly 
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invertible. Its inverse function T(v) is differentiable, hence so is W(v) = 
w(T(v)), and these variable satisfy the nonlinear differential equation 
The explosion time (resp. the value of the adaptation variable at this time) 
is the limit ofT{y) (resp. W{v)) when v —¥ oo. 

Therefore the description given by equations ^ is valid in particular during 
the explosion of the membrane potential variable. Both the time of the spike 
and the value of the adaptation variable at this time can be accurately and easily 
simulated by solving the smooth ordinary differential equations (]6]), that governs 
the trajectory W(v) of the adaptation variable in the phase plane and T(v) the 
inverse function of t h-> v(t). The description of the trajectories in the phase plane 
allows simulating well-behaved equations during this phase where the Euler scheme 
was particularly sensitive. This description is the cornerstone of the algorithm we 
develop. 

The discrete integration scheme makes use of this convenient description of 
the trajectory in the phase plane in regions of the phase plane where the function 
t i — y v(t) is smoothly invertible, and in the rest of the phase plane close to the 
f-nullclines, simulates the solutions of the equations Q with a fixed time step 
Euler scheme. 

We propose two algorithms based on this idea. The first one is a fixed inte- 
gration step method and the second one a fixed precision method with adaptive 
integration steps. 

3.1 Fixed integration-step method 

Let us assume that we have defined M > a real parameter, and let dt > and 
dv > be respectively a time step and a space step that will be used to numerically 
compute the solution of equation Q and (pi. 

Starting from an initial condition (vq, wq) at time to, we build iteratively a 
sequence (t n ,v n ,w n ) approximating the orbit (t,v(t),w(t)) solution of the contin- 
uous ordinary differential equations Q. Given (t n ,v n ,w n ) the numerical solution 
computed at iteration n, we compute t n+ %, v n+ i and u>„+i as follows: 
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(i) . If \F(v n ) — w n + I(t n )\ < M, then we update the values of the triplet as: 

v n+1 =v n + dt (F(v n ) - w n + I{t n )) 
< w n+1 = w n + dta(bv Tl - w n ) 

tn+l — t n + dt 

corresponding to the integration of equation ([I]) with a constant time step 
dt. This phase of the integration algorithm is referred as the time-integration 
phase. 

(ii) . if \F(v n ) — w n + I{t n )\ > M, then we update the triplet as: 



< 





= V n + dv 


W n+ l 


= w n + dv 


tn+1 


= t n + dv - 



a (bvn—w n ) 
F(v n )-w„+I(t n ) 



F(v n )-w n +I(t n ) 



which corresponds to the integration of (J6J) with constant w-step dv. This 
phase of the integration algorithm is referred as the phase-plane integration 
phase. 

Intuitively, this method consists exactly as announced in solving ([T]) in the re- 
gions of the phase plane where the vector field is of small amplitude and solving ^ 
in the regions of phase plane where the vector field of the dynamical system given 
by ([T| is of large amplitude and yield divergence in finite time of the solutions. 

Let us now analyse this algorithm and control the precision of the solution 
computed. First of all, we analyse the solutions of the algorithm in spiking regions 
when no threshold is taken into account. 

Lemma 3. For any initial condition (v ,Wq) in the spiking zone, the variable v n 
blows up in finite time. 

Proof. Let (t>o,wo) £ Z*, we have for any n, (v n ,w n ) G Z*. After a finite number 
of iterations K, we will have for all n > K F(v n ) —w n + I(t n ) > M. Indeed, while 
(v n ,w n ) E Z*, lemma [I] implies that v n is unboundedly increasing and moreover 
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in that zone we have 



F(v n ) -w n + I{nr) > F(v n ) -bv n + P. 

Since v n is unboundedly increasing when simulated with Euler scheme, so is 
F(v n ) — bv n + I (tit) and therefore there exists K G N such that M < F(vk) — 
bv K + I* < F(v K ) -w K + I(Kt). 

The next step vk+i, wk+i is hence computed using the simulation in the phase 
plane algorithm. Using the same argument as in the proof of lemma [TJ we obtain 
that the sequence (v n , w n ) stays in a zone where F(v) — w + I > M and w < bv, 
and therefore for all n > K, the simulation algorithm used is the phase plane 
algorithm. 

We therefore obtain that v n tends to infinity since t> n+1 = v n + dv. We now 
show that the sequence t n converges to a finite limit. Indeed, we have for all 
n > K: 



n 

<t K + dv 



, v _ i F(v K + (n-K) dv) -b(v K + (n- K) dv) + I* ' 

The series converges when n — > oo because of the hypothesis that F grows faster 
than v 1+s for some 5 > 0, by a simple sum/integral comparison. □ 

In order to compute the precision of the algorithm in the evaluation of the 
spike time and on the value of the adaptation variable at the time of the spike, 
we now consider a thresholded version of the simulation algorithm consisting as 
usual in introducing a cutoff 9. The simulation algorithm is run as long as v n < 9. 
When v n exceeds 9, we instantaneously reset the membrane potential and update 
the adaptation variable by defining: 



v n +i = c 
w n+1 =w n + d 

tn+l t n 
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Figure 2: Partition of the phase space. Pink zone: phase-plane integration and 
white zone: time integration. Orange dotted box: initial conditions. 



We are now in position to to propose a suitable choice of (M, dt, dv) ensuring 
a fixed precision e on the adaptation variable at the time of the cutoff v = 9. An 
optimal choice of (M, dt, dv) would moreover minimizes the number of operations 
necessary to reach the desired precision e, but this optimal choice is quite diffi- 
cult to define. We nevertheless define an even more efficient method in section 



3.2| that involve adaptive integration steps ensuring a fixed absolute error e. To 
this purpose, let us assume that we consider a bounded set of initial conditions 
[i>;l,i>m] x [w m ,WM} (orange box of figure Fig. [2J, and let us consider the regions 
of the phase plane reached by the orbits of the system starting from these initial 
conditions. 

Lemma 4. For any initial condition (v , w ) G [v m , %] x [w m , Wm\, the thresholded 
trajectories (v(t),w(t)) are contained in a compact set [vi,9] x [wd,w u ]. Moreover, 
if F grows faster than v 2+s for some 5 > 0, then v\, and w u are independent of 
9, and the time integration phase takes place in a compact set independent of 9. 

The proof of this quite intuitive lemma is provided in appendix [Bj and is 
illustrated in Figure [2] the zones of the phase plane depicted in gray, and the rest 
of the phase plane not plotted in the figure constitutes a zones of the phase plane 
the orbit starting from any initial condition in the subset defined will never reach. 
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Now that we restricted possible initial conditions of the system, we know from 
lemma [4] that the values of (v,w) on the orbits are contained in a bounded set. 
For the quadratic model, these values depend on the cutoff value 6 since the 
adaptation variable at the times of the spike diverges, and for models satisfying 
assumption 



(A2) , it can be defined independently of 6. On this bounded set, 



the functions F(v) — w + I and a(bv — w) are both bounded, differentiable with 
bounded derivatives, which allows definition of time steps dt and phase-space step 
dv that uniformly ensure a given precision of the numerical algorithm. This fact 
is discussed in appendix [Bj and we now turn to discuss a more efficient method 
based on the algorithm proposed but involving adaptive integration steps. 

This algorithm can be adapted to most classical simulation tools used in com- 



putational neuroscience (e.g. BRIAN (Goodman and Brette 2008) or NEST (Dies- 



mann and Gewaltig, [2002)) and that make use of fixed time step methods, by 



including a condition to switch either to integrating (i) or (ii) 



3.2 Adaptive integration steps algorithm 

Let us consider a constant M > fixed. In contrast to the fixed integration step 
methods introduced in the previous section, we define e the absolute precision we 
wants to achieve all along the simulation of the orbits. We also define DT and DV 
the maximal values of time and phase-space integration steps we want to take into 
account. An efficient method minimizing the computation operations consists 
in defining at each operation the optimal integration steps ensuring a precision 
bounded by e at each step of the integration algorithm. The absolute precision 
is known to be bounded by the second derivative of the vector field multiplied by 
the integration step. We will therefore define the integration step as a function 
of the value of the modulus of the derivative of the vector field in each case. Let 
(to,vo,wo) be an initial condition of the dynamical system (JT1) with cutoff and 
reset, and assume that we iteratively computed the solution up to rank n. We 
therefore dispose of a point (t n ,v n ,w n ). We compute the next point of the orbit 
(Wii v n+1 , w n+1 ) as follows: 

(i). If \F(v n ) — w n + I(t n )\ < M, we simulate the solutions of equation Q. The 
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time step chosen depends on the value of the second derivative of v(t) and 
w(t), that read: 

v "(t) = F'(v)(F(v) -w + I)-a(bv-w) + I'{t) 
w "{t) = a (b (F(v) - w + I))-a(bv-w) 

where x' denotes the derivative of x(t) with respect to time. The larger time 
step dt ensuring a precision of order e at the point (v, w) therefore reads: 

dt(v, w) 



max(|t»"(t)|, \w"(t)\)' 



In order to keep enough points in regions where the vector field varies slowly, 
we define r n = min(dt(v,w), DT) and we define: 



v n+1 =v n + r n (F(v n ) -w n + I(t n )) 
w n+1 = w n + r n a(bv n - w n ) 

tn+l tn + T n 



(ii). If |F(t> n ) — w n + I(t n )\ < M, we advance the simulation by computing the 
next point through equation ([!]). The integration step chosen depends on the 
value of the second derivative of T(v ) and W(v ) with respect to v (denoted 
with a double dot), that reads: 



W(v) 



T(v) 



a b 



a{bv-W) aF'{v){bv -W) 



F(v) -W + I (F(v) — W + I) 2 (F(v) -W + I) 2 
a(bv-W)(a(bv-W) + I'(T)) 

(F(v) - W + I) 3 
_F'(v)(F(v) -W + I)-a(bv-W) + I'(T) 

(F(v) -W + I) 3 



With these expressions, we are in position to choose an optimal integration 
step depending on the values of the variables (v, w) ensuring the desired 
precision: 

dv(v,w) = r. — r. . 

V ; max(\W(v)\,\T(v)\) 



24 



Here again, in order to keep enough points in regions where the vector field 
varies slowly, we define dv n = mm(dv (v , w), DV) and we define: 



v n +i =v n + dv n 



„, , - 7/1 4- fiii a(bv n -w n ) 

w n+1 -w n + av n F(Vn) _ Wn+I(tn) 



tn+i — t n + dv r , 



i 



F(v n )-w n +I(t n ) 



The integration steps chosen at each pace of the algorithm ensures that the 
precision is always bounded by e. During phase (i), it therefore ensures that we 
have 

max{\v(t n ) - v n \, \w(t n ) - w n \) < e 
and during phase (ii), that 

max(\W(v n ) - w n \, \T(v n ) - t n \) < e. 

This algorithm therefore ensures an overall precision bounded by a constant C 
multiplied by e, where the constant C is proportional to the maximal derivatives 
of v i — y W(v) and v i— > T(v), both of which are uniformly bounded on Z*. 

Note that the choice of the value of M impacts the number of operations 
performed. Choosing a large value for M implies that the time integration phase 
is emphasized, whereas a small value of M emphasizes the phase space integration. 
Note also that the remark on the explosion of v in the algorithm is still valid, and 
that on a spiking trajectory, the integration step in phase (ii) of the algorithm is 
soon chosen to be DV. 



3.3 Numerical Errors accumulation during a spike train 
simulation 

Similarly to the case treated in section |2j we have shown that the evaluation of the 
spike time and of the value of the adaptation variable at this time yielded errors, 
but in contrast with fixed time-step method, the absolute precision is controlled in 
our algorithm. Indeed, let us assume that (v(t),w(t)) is a spiking orbit, and that 
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for t > £*, the orbit is contained in Z* and moreover F(v(t)) — w(t) + I(t) > M. 
Let us denote by X* = (t*,v*,w*) a point of the orbit in this region of the phase 
space and X** = (t**, v**, w**) another such initial condition. We are interested 
in the divergence of the trajectories starting from X* and X** and solutions of 
equations (|6]), denoted respectively (Ti(v), W\(v)) and (T 2 (v), W^{v)). Integrating 
formally the equation and subtracting the expressions obtained, we get, denoting 
Gi = (F(v) - W x {v) + /(Tx^))) and G 2 = (F(v) - W 2 (v) + I(T 2 (v))) 

T x {v) - T 2 (v) = (t - t**) + £7 ^du + /;„ a^-a^du 
W x {v) - W 2 (v) = (w* - w**) + £7 a(bU aZ U)) du + X» aibU G^ U)) - a ' bU aZ U)) du 

These equations are contracting and will result in an exponential decrease of the 
error made. Let us demonstrate this fact by first considering the case where the 
input current I is constant (the general case can be done along the same lines). 
Making more explicit the above expressions, we get: 



Ti{v)-T 2 (v) = {t*-t**) + 



1 du+ rm^md« 



Gi(u)G 2 (u) 



and therefore the error made on the evaluation of the spike time is governed by 
the error made on the evaluation of the adaptation variable, which reads: 



W x (y)-W 2 (y) = (w*-w* 



V* 



a(bu-Wi{u)) 



du- 



v a{bu-{F(u) + I)) 
G\{u)G 2 {u) 



(Wi-W 2 



Let K > fixed. There exists L > v** such that for all u > L we have (bu— (F(u)+ 
I)) < 0. Therefore, as v increases, the equation on the distance between W\ and 
W 2 becomes contracting and the distance between the two adaptation variable 
trajectories is reduced. This result extends to the time variables, and hence in 
contrast to the system ([T]), the error is here reduced instead of exponentially 
increased, and the larger the cutoff 6 is, the smaller the distance between the two 
computed trajectories^] 

^The complete proof of this result involves Gronwall theorem, and can be related to the 



contraction theory (see e.g Slotine and Li (1991); Lohmiller and Slotine (19981) 
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After the first spike is fired, the numerically computed spike time t** and 
the value of the adaptation variable at this time w** differ from the actual spike 
time t* and value of the adaptation variable w*, and the difference is bounded 
by the precision e chosen. We have noticed that the equation ([TJ was expanding 
along spiking trajectories, i.e. that the distance between two orbits starting from 
different initial conditions diverge close to a spike emission. However, in the case 
of the algorithm we propose here, instead of simulating equations close to the 
spike emission, we rather simulate equations We showed that these equations, 
contrarily to the original dynamical system, are contracting when v is large enough: 
two orbits starting from different initial conditions converge exponentially the one 
towards the other. Therefore, instead of increasing the errors, the simulation of 
will reduce the error made on the evaluation of the spike time and on the 
adaptation variable at this time. Though the discontinuity might produce large 
errors at the reset points, the eventually computed spike time and adaptation 
variable are reliable because of the contracting properties of the flow of equations 




3.4 More refined integration procedures 

In the algorithm we proposed in this section, we made use of a Euler scheme, either 
with fixed or adaptive integration steps. The exact same procedure can be per- 
formed with more refined numerical integration procedures, such as Runge-Kutta 
methods, Adams-Moulton or Adams-Bashforth multistep methods, or implicit 
schemes, in each phase of the algorithm, namely the phase (i) of time integra- 
tion and the phase (ii) of phase-space integration. These methods will therefore 
present both the advantages of finer integration methods and of our integration 
scheme that chooses between the two description of the orbits and (|6| the 
one that minimizes numerical errors, and also the fact that the errors made on the 
evaluation of the spike time and the adaptation variable at this time exponentially 
decreases in the integration phase (ii). These therefore provide a set of efficient 
simulation methods. 
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4 Discussion 



In this paper we studied the error produced by fixed time step integration methods 
for nonlinear bidimensional neuron models and showed that the error produced, 
even if these are of order 1 in the time step chosen, are proportional to a function 
of the parameters and of the cutoff chosen to identify the spikes, and this function 
is unbounded as the cutoff is increased. This fact implies that in order to achieve 
a fixed precision on the global algorithm and accurately simulate the explosion 
time, one needs to use a large enough threshold value 9 and a very small time step 
in order to accurately estimate the spike time and the value of the adaptation 
variable at this time. We saw that the main difficulty arises close from the spike 
time, since the explosion of the membrane potential variable produces large errors 
and the differential of the adaptation variable also diverges at this time. The 
small time step needed to accurately evaluate the crucial values of the spike time 
and the adaptation variable at these time will result in an excessive precision in 
regions where the vector field is smooth and an increased computational time. 

In order to overcome this difficulty, we proposed an alternative algorithm that 
is based on the numerical computation of the orbits in the phase plane in regions 
where the modulus of the vector field is large enough. This algorithm, similar 
to classical variable integration step methods, allows to choose optimal steps in 
order achieve the desired precision in a minimal number of operations, which can 
be uniformly bounded as shown in section [3j Therefore this method provides us 
with a stable method that allows considering large thresholds and that precisely 
evaluates the spike times and the value of the adaptation variable at these times 
in a reduced number of operations. 

In this section we specifically focus on the impact of numerical errors in the 
evaluated trajectories and on the computational efficiency of the algorithms. 

4.1 Quantitative imprecision yield dramatic qualitative er- 
rors 

We now discuss the importance of numerical errors on the simulations from a 
computational neuroscience viewpoint. The models we addressed in this paper 
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are particularly used to simulate spike times in order to identify different spike 



patterns, such as regular spiking, bursting and chaotic spiking (see e.g. (Izhike- 



vich ( 2003 2007 Brctte and Gerstner, 2005 Touboul 2008a Touboul and Brette 



2009)). We showed that the spike times and the adaptation variable at these times 



can be substantially modified when using fixed time step methods, and this effect 
can produce qualitative distinctions between firing patterns. Indeed, it was shown 



in (Touboul and Brette, 2009) that the spike pattern depended sensitively on the 
value of the adaptation variable at the times of the spikes and that these spike 
patterns undergoes bifurcations. Since we have seen that fixed time-step Euler 
scheme produces a systematic negative error (i.e. the evaluated spike time was 
strictly inferior its actual value as a solution of the continuous differential equa- 
tions), this error can make the system cross bifurcation points and change spike 
patterns. 

We illustrate this fact for instance on the quadratic model with original pa- 



rameters provided in (Izhikevich, 2003): in this case, F{x) = 0.04a; + 5x + 140 
and c = —59.9, parameters that were chosen to fit intracellular recordings. By 
choosing the parameters a = 0.02, b = 0.19, d = 1.15, 9 — 30 and a constant input 
/ = 7.6, the system produces bursts with two spikes per bursts. We simulate for 
T = to 1000 in order to record a long enough spike train, and plot the sequence of 
values of the adaptation variable at the time of the spikes. These values allow dis- 
criminating between bursts (that produce periodic sequences of reset values with 
period two) and regular spiking (that correspond to a fixed point of the sequence). 
Results of the simulation are plotted in Figure [3j We observe that for dt < 0.025 
the bursting nature of the trajectory is recovered, but for any time step greater 
than this value, the bursting nature of the spike pattern produced is not clearly 
identifiable. Moreover, we observe that for dt > 0.05, the spike patterns produce 
corresponds to a regular spiking behavior: the system crossed the period-doubling 
bifurcation because of numerical errors, and this is not a purely visual effect: the 
sequence of reset values is unimodal (as shown in Figure [3^ A)). Therefore, on our 
simulation interval, one needs to perform at least 30 000 operations to recover the 
nature of the desired spike train, and at his resolution the burst is quite perturbed 
by numerical imprecision. 
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Figure 3: Simulation of the quadratic model with parameters given in the text. 
(A) represents the value of the adaptation variable at the spike times as a function 
of a fixed time-step dt for an Euler scheme, (left) We observe that the bursting 
nature of the trajectory is lost and one recovers a regular spiking behavior if the 
time step is not small enough, (right) represents the histogram of these values 
for dt = 0.01 (green) and dt = 0.1 (red). (B) and (C) correspond to sample 
trajectories computed for dt = 0.01 and dt = 0.1 respectively (stopped at t = 
250 for legibility). (D) represents the sequence of reset values for the variable 
integration step algorithm with precision set at e = 0.1 (similar to (A) but for 
a fixed integration scheme and where the abcissa corresponds to the number of 
spike fired before this value is obtained) and (E) a sample trajectory simulated by 
this algorithm. 
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For the same parameters, we set a precision e = 0.01 and compute the trajec- 
tory using our fixed precision algorithm. The algorithm performs 2000 operations 
and achieves a great precision on the value of the adaptation variable at the spike 
time, that is not comparable to what is obtained by the fixed time-step Euler 
method. Actually the precision of the algorithm is smaller than the value of e 
chosen because of the control put on the value of the integration steps in order for 
these not to exceed certain values in regions where the vector field varies slowly. 
For such precision, one would necessitate to use dt = 0.01 corresponding actu- 
ally to 100 000 operations. Therefore, the computational gain of using the new 
algorithm instead of a fixed time step Euler method ranges from 30 to 50. 

We now compare our algorithm to different classical algorithms with adaptive 
time-steps encoded in Matlab. All classical algorithms are tested, and both sim- 
ulated trajectory and execution time are compared to our algorithm. For small 
time intervals, all methods fairly well approximate the solution, but after some 
time, the spikes are shifted and though all methods still conserve the bursting 
nature of the trajectory, all the solutions are increasingly shifted. We set an abso- 
lute precision of 0.01, and compare the computation times, summarized in table [TJ 
and the trajectories. We observe that the algorithm proposed in the manuscript 
is way faster than any other method, which was expected since it is adapted to 
the very nature of the problem and is not as generic as Matlab's routines. We 
observe that all the simulations increasingly large delays in the spikes compared 
to our simulation algorithm (see Figure [1]), which is probably linked with the fact, 
evidenced here in fixed time-step Euler integration scheme, that estimated spike 
times with fixed time step integration schemes or adaptive time step algorithms 
present a systematic positive delay in the estimation of the spike time, and also 
that the errors accumulate all along the emission of a spike train. However, an 
adaptive method seems to stand out: the ode 113 routine, implementing a vari- 
able order Adams-Bashforth-Moulton PECE solver, a linear multistep integration 



method (see e.g. (Hairer et al. , 1993)) provides a very good match with our simu- 
lation algorithm, with reasonable computational time. None of the other Matlab's 
algorithms reliably simulate the system, even if the precision is set to e = 0.01, 
on the interval [0, 1000]. This will also be the case of the fixed time-step Euler 
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Algorithm 


Paper's 


ode23 


ode45 


odell3 


odel5s 


ode23s 


ode23t 


ode23tb 


Exec, time 


0.0053 


0.7478 


0.6730 


0.6567 


1.993 


2.301 


1.7538 


1.7756 



Table 1: Compares the execution time of the paper's algorithm with different 
algorithms proposed in Matlab. 
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600 650 700 750 800 850 900 950 1000 



(a) All Algorithms 



(b) odell3 



Figure 4: Comparison of the solutions computed with different integration meth- 
ods implemented in Matlab. The variable order Adams-Bashforth-Moulton PECE 
solver (ode 113) provides a very good match with the algorithm proposed in this 
paper. 



scheme: whatever the time step dt is, small variations will progressively shift the 
spikes. Besides this imprecision, all algorithms present moreover longer execution 
time evaluated on a Macbook with processor 2 GHz Intel Core 2 Duo and mem- 
ory 2GB 1067 MHz DDR3. However, though these differences, it is important to 
note that all the algorithms reproduce the bursting property and provide a good 
match with the sequence of reset values expected. This is not the case of the 
Euler method for time step dt > 0.03. This conclusion extends to most adaptive 



time-step integration methods, such as the one implemented in NEST (Diesmann 



and Gewaltig, 2002). 



4.2 Number of operations needed to achieve fixed preci- 
sion for fixed time-step Euler scheme 

In this section we compare the computational efficiency of the two proposed algo- 
rithms. In Figure [TJ we plotted the error made on the evaluation of the adaptation 
variable as a function of the cutoff 9 chosen. Let us denote by E(9) this function. 
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The error is then proportional to E{9) dt and therefore these curves are directly 
related to the time step needed to achieve a given precision. Indeed, if one wants 
to get an absolute precision bounded by e on the adaptation variable, he needs 
to choose dt < e/\E{9)\ which implies that the number of operations needed is 
proportional to \E(6)\/e, which diverges as the cutoff is increased. The curves on 
the error |A(i>)| and are therefore proportional to the number of operations 

needed to achieved fixed precision. 



4.3 Further Topics 

In models satisfying assumption (A2)| such as the quartic and the exponential 
adaptive models, the adaptation variable at the times of the spike converges, and 
therefore it can be very convenient to use large values of the cutoff in order to get 
rid of the artificial strictly voltage threshold similar to classical integrate-and-fire 



models (Lapicque, 1907 Stein, 1967) by a more realistic spike initiation (see e.g. 



(Latham et al. 2000 Fourcaud-Trocme et al. , 2003)) and fully taking advantage 
of this property of nonlinear integrate-and-fire neurons. The algorithm proposed 
here allows considering larger values of cutoff 9 without a prohibitive increasing 
the computational time, which allows more accurately simulating the blowing-up 
system. Another interest of the proposed numerical method is the flexibility in 
terms of spike times: in contrast to all fixed time-step methods where the spike 
times are aligned on a grid, our algorithm can have arbitrary spike times, which 
is an important property from the biological viewpoint and for the computation 
of large networks or over large period of times. 

The algorithm we proposed here can easily be used for simulations of networks 
composed of neurons of the type studied here. This new method will therefore 
be very useful for the simulations of very large networks using such models by 
providing a more precise and computationally more efficient method to simulate 
the solutions of these equations. 
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A Some Theoretical Results 



In this appendix we review and prove some properties of the dynamical system 
Q of interest in the main text. 

We first recall that this dynamical system either has zero, one or two fixed 
points depending on the values of the parameters. More precisely, let us denote 
by m(b) the unique minimun of the convex function v \- > F(v) — bv, and by 
v*(b) the unique point where this minimum is reached (i.e. the unique solution 
of F'(v*(b)) = b, which exists under assumption (Al)[ ). It is proved in (Touboul 



2008aP (see figure |5J) that: 

• if / > —m(b) then the system has no fixed point; 

• if / = —771 (6) then the system has a unique fixed point, (v*(b), w*(b)), which 
is non-hyperbolic. It is unstable if b > a. 

• if / < —m(b) then the dynamical system has two fixed points b), w_ = 
bv-(I,b)) and (v+(I, b), w + = bv + (I, b)) such that 

v-(I,b) < v*(b) < v+(I,b). 

The fixed point v+(I, b) is a saddle fixed point, and the stability of the fixed 
point f-(J, b) depends on I and on the sign of {b — a): 

(i) . If b < a then the fixed point v_(I, b) is attractive. 

(ii) . If b > a, there is a unique smooth curve I s (a,b) defined by the im- 

plicit equation F'(v-(I s (a, b),b)) = a. This curve is given by the equa- 
tion I s (a,b) = bv*(a) — F(v*(a)) where v*(a) is the unique solution of 
F'(v*(a)) = a. 

(a) If / < I s (a, b) the fixed point is attractive. 

(b) If / > I s (a, b) the fixed point is repulsive. 



In (Touboul and Brette, 2009), the authors thoroughly define a Markov parti- 
tion of the dynamics in each case, that involve a spiking zone which is stable under 
the flow and where any orbit with initial condition in the zone blows up in finite 

34 



(A) 

t (B) 




Figure 5: (A) Fixed points and stability in the parameter plane (I,b). The pa- 
rameter zone PI corresponds to a case with no fixed points, and the related zone 
Z is displayed in (B). In the parameter zones P2 and P3, the Z zone is displayed 
in(C). 

time (and fires a spike). The precise description of this zone is quite complex and 
involves the description of the stable manifold of the saddle fixed point v + (I,b), 
which does not have a simple analytical expression. We are interested here in 
defining a sub-zone of the spiking zone stable under the flow of the differential 
equation that allows a simple computation of the spike time and of the adaptation 
variable at this time. In definition [TJ we introduced a subset of the phase plane 
Z*. This subsets of the phase plane present the interest to satisfy the following 
properties 

Theorem 5. Assume that the input current I(t) depends on time, and moreover 
that t i — y I(t) is lower bounded (at least on the time interval considered (to,Ti)), 
i.e. I(t) > I* for any t G (t q ,Ti). Then we have: 

(i). Z(I*,b) =: Z* is stable under the flow of the equation 

(ii). Let (v ,w ) G Z* , and denote (v(t),w(t)) the solution of equations ([T| hav- 
ing initial condition (vq,wq) at time to G [to,ti]. Then t i— > v(t) is strictly 
increasing, therefore invertible. Its inverse function T(v) is differentiable, 
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hence so is W(v) = w(T(v)), and these variable satisfy the nonlinear differ- 
ential equation: 

dT 1 



dv F(v) -W + J(T) 

(7) 

dW a(bv-W) 

~ F(v) -W + I(T) 

The explosion time (resp. the value of the adaptation variable at this time) 
is the limit ofT(v) (resp. W(v)) when v — > oo. 

Proof, (i). The stability of Z* = Z(I*,b) under the flow of the equation stems 
from the fact that the vector field points inwards in Z* on its boundary (see 
figure [5]). Let us first address the case I* > —m(b). In that case, we have 
(by definition of m(b)) F(v) -bv + I(t) > for all v G R and all t G (t ,Ti). 
Z* is defined as {(v,w) ; T(v,w) < 0} where T(v,w) = w — bv. We use a 
proof by contradiction. Let us assume that the zone Z* is not stable under 
the flow, and that a solution (v(t),w(t)) with initial condition (vq,wo) G Z* 
at to G (to,Ti) exits Z* at time t\. This time t± is defined as the first exit 
time of Z*, i.e. 

t 1 =in£{t>t ; (v(t),w(t))iZ*}. 

Because of the continuity of the boundary and of the solutions, necessarily 
we have w(t\) = bv(t\). Moreover, at this point we have: 



dT(v(t),w(t)) 



dt 



= a (bv(t) - w(t)) - b (F(v(t)) - w(t) + I(t))\ t=ti 

t=ti 

= -b (FivfaV-bvitJ+Ifa)) 

< -biFivihV-bvi^ + T) (8) 



This last expression is strictly negative because of the fact that I* > —m(b). 
Therefore, since 7 : t 1— > T(v(t),w(t)) is a differentiable function of time 
such that 7(ti) = and ^1*=^ < 0, there exists 77 > such that 7(2) < 
on [ti,ti + 77), hence (v(t),w(t)) G Z* on + rf\ which contradicts the 

definition of time t\ as the exit time from Z*. 

In the case where /* < —m(b), the zone Z* is defined by the additional 
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condition that v > v+(I*,b). On this zone, v h-> F(v) — bv is non-negative 
and increasing because of its convexity property. Therefore, on Z*, we have 
(F(v(ti)) — bv(ti) + I*) > 0. Let us consider (v(t),w(t)) a solution with 
initial condition (v ,w ) G Z*. We have v > v + (I*,b). While (v(t),w(t)) is 
in Z*, v(t) is strictly increasing, and hence the solution will never cross the 
boundary v = v + (I*, b), and can only exit Z* through the boundary w = bv. 
On this boundary, the same argument as in the case where J* > —m(b) 
applies, which ends the proof. 

Note that when the current I(t) is constant, this result directly comes from 
the property that the vector field points inwards the zone Z* (and inwards 
any Z(V) for V > v+(I,b), with the convention v + (I,b) = — oo if I > 



-m(b)), making of the zone Z* a trapping region as referred in (Strogatz 



1994, Chap. 7.3.) (a proof of this property can be found for instance in 



flViterbol|2005) ) 



(n). Let (v ,Wq) G Z*. Then necessarily (v ,w ) G Z(V) for some V such that 
v+(I, b) < V < v . In Z(V), we have F(v) — w + I > F(v) - bv + J > 0, 
therefore t i— > v(t) is a strictly increasing function. Moreover, since the 
derivative of v(t) never vanishes, t i-> v(t) is invertible with differentiable 
inverse function T(v) by application of the inverse function theorem. 

From the stability property of Z*, for all t >to we have w < bv and hence 
^1 = F(v) -w + I{t) > F(v) -bv + P 

{Jib 



whose solution blows up in finite time under assumption (Al) , by the virtue 



of Gronwall's lemma, see (Touboul, 2009). 



We denoted by T(v) the inverse of the function t h-> v(t) and noticed that 
it was differentiable, and therefore W(v) the composed application W(v) = 
w(T(v)) is also differentiable. It is clear that the explosion time corresponds 
to the limit of T(v) when v — > oo and the value of the adaptation variable 
at the explosion time the limit of W(v) when v — >■ oo. We have by definition 
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v(T(v)) = v, and differentiating by v we get: 



dv(T(v)) dv(t) 



dv 



dt 



dT(v) 



t=T(v) 



dv 



(F(v(T(v)))-w(T(v)) + I(T(v))) 



dT(v) 
dv 



and it is also equal to one since v(T(v)) = v. Therefore we have, since 
F{y) — w + I(t) never vanishes on Z(V): 

dT(v) 1 



dv F(v)-W(v) + I(T) 



and using the differentiation formula of a composed function, 



dW 



a(bv-W) 



dv F(v)-W + I(T)' 



□ 



Another property we can deduce from (Touboul, 2009) is that the adaptation 



variable in a left neighborhood of the explosion time is always negligible compared 
to v. This remark implies that the differential of w(t) tends to infinity at the 
explosion times of v. However, theorem [2] shows that in the region where the 
neuron elicits a spike, the time of the spike and the adaptation variable satisfy a 
well-behaved differential equation ([7|), fact on which our simulation algorithm of 
section [3] is grounded. 



B Trajectories starting from a given set of initial 
conditions are bounded 

In this appendix we demonstrate lemma |4j which is recalled here: 

Lemma 6. For any initial condition (vq, wq) G [v m , %] x [w m , wm], the thresholded 
trajectories (v(t),w(t)) are contained in a compact set [vi,9] x [wd,w u ]. Moreover, 
if F grows faster than v 2+s for some 5 > 0, then vi, Wd and w u are independent 
of 6, and the time integration phase takes place in a compact set independent of 
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6 (see Figure^ the zones of the phase plane depicted in gray, and the rest of 
the phase plane not plotted in the figure constitutes a zones of the phase plane the 
orbit starting from any initial condition in the subset defined will never reach). 

Proof. Let (t>o,u>o) £ [v m ,VM] x [w m ,WM] and consider the orbit of the dynamical 
system ([I]) with (v , wq) an initial condition in this region. We start by considering 
I(t) constant. The trajectories for 1(f) non constant will be bounded by the values 
the orbit can reach for the vector associated with the minimal and maximal values 
of the I(t) because of the monotony of the vector field and because of Gronwall's 
inequality. 

First of all, the case where F(v) — bv + / > for all v corresponding to the 
case where the nullclines never intersect is simple to treat. Indeed, in that case 
the spiking zone is simply {(v,w);w < bv} and it is a trapping zone, i.e. any 
trajectory with initial condition in the spiking zone will never escape from it. For 
any initial condition in this spiking zone (zone I in Figure 

gA)),w < w(t) <bv(t) 

and v < v(t) < 6 and hence w(t) G [w m ,max(b9, wm)] and v G [v m , 9}. If (v ,w ) 
in zone II of Figure gA), i.e. if bvo < wq < F(vq) + I, the variable v(t) is 
increasing all along the trajectory and therefore 6 > v (t) > vq > v m . The variable 
w(t) is initially decreasing (therefore bounded by wq < wm) before reaching the w- 
nullcline and then increases in the spiking zone and where trajectories are bounded 
by b9. Hence bv m < bvo < w(t) < max(wjw, bd). Eventually in zone III of 
Figure gA), i.e. w > F(v ) + I, both the value of v(t) and of w(t) are initially 
decreasing and will reach in finite time the v-nullcline where it enters zone II. 
Before reaching this zone, the minimal value reached by the variable v(t) is greater 
than V\ = min(F~ 1 (wo — I)) and smaller than Vq < Vm, and the value of w(t) is 
smaller than wq < wm and greater than min(F(f ) + 1) which is reached at a point 
denoted v* . From these initial conditions in zone II the same analysis as previously 
done extends the boundaries to 9 > v(t) > min(t> m ,t>i) and b mm(vi,v m ,v*) < 
w(t) < max(w M ,6^). 

After resetting, the initial condition is on the line c with values of w in the 
interval [mm(vi,v m ,v*) + d, max(wM,b6) + d] and a similar analysis provides us 



39 



with the boundaries 



(v(t),w(t)) E [min(t) m , n*)] x [b min(t;*, v m , v*, c), max(w A/ , 

with t>* = min(i 71_1 (max(iL'Af , b9) + d — I)). 

Moreover, in the case where F grows faster than v 2+6 for some 5 > 0, it was 



proved in (Touboul and Brette, 2009) by a thorough analysis of the spike and 
reset process that the maximal value reached by the variable w in the spiking 
zone is bounded by a value <3>* independently of the cutoff, making the boundaries 
provided independent of 9. 

In the case where the nullclines do meet, we perform a similar analysis. First of 
all, it is clear that the value of w(t) before resetting is bounded by max(wM, b9). 
After the reset, w is therefore upperbounded by w u = max(w M , b9) + d. The value 
of v is cut at 9 and hence any thresholded trajectory have v(t) < theta. We assume 
that F(c) — w u + I* < 0. At this point we therefore have v decreasing, and the 
minimal value of v after a reset corresponds to this value and is lowerbounded by 
the solution of F(v) - w u + F = i.e. v > mm(F~ 1 (w u - J*)). If w u < F(c) +F, 
then the value of v on any reseted trajectory is greater than c. Therefore, we 
conclude that v > v\ = min(c, v m , min(F -1 (u7 u — /*))). Eventually, this value 
provides us with a minimal value of w along the trajectories. Indeed, the minimal 
value of v after reset implies that any reseted trajectory will have w(t) > bv\. If 
w is such that F(v ) — w Q + I* < 0, then in this region 

First of all, because of the structure of the flow and the reset condition, the 
maximal value reached on the reset line v — c is upperbounded by b 9 + d. This 
bound can be made independent of 9 if the function F grows faster than v 2+6 
for some 5 > 0. In that case, there exists a value denoted $* + d, as defined in 



(Touboul and Brette, 2009), and corresponding to the value where the adaptation 
variable will be reset when starting just below w* = F(c) + I*. This value does 
not depends on theta and can be evaluated previously before simulations, by 
simulating equation ([7]) with initial condition (c, wq) with wo slightly below w* = 
m(b) the minimum of the function F(v) — bv + 1. Note that if 9 is not chosen 
very large, or if one does not want to vary 9, this value $(iy*) can be for instance 
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replaced by the an upper-bound b9 + d (choice depicted in the figure). Therefore 
the values of the adaptation variable at the time of the spike have the global upper- 
bound w u = max(%, <3>(u>*) + d) or w u = max(%, b6 + d). This maximal value 
defines a left boundary of the values of v. Indeed, denoting v% = mm{v, F(v) — 
w u + I* = 0} (there are two possible values of W\ because of the convexity of F), 
we have v > vi := min(f!,w m ). This value in turn allows defining the lower w 
value trajectories can have. As soon as the w-nullcline is crossed, w increases. 
Therefore, w is larger than wi := mm(bvi,w m ). It has a typical structure as 
depicted in figure Fig. |2} The set where the algorithm (i) is used is the white 
space in figure [2j and is strictly included in [vi, v+(I* — e, b)] x [wi,w u ]. □ 

Now that we restricted possible initial conditions of the system, we know from 
lemma [4] that the values of (v,w) on the orbits are contained in a bounded set. 
For the quadratic model, these values depend on the cutoff value 6 since the 
adaptation variable at the times of the spike diverges, and for models satisfying 



assumption (A2) it can be defined independently of 6. On this bounded set, 
the functions F(v) — w + I and a(bv — w) are both bounded, differentiable with 
bounded derivatives, which allows definition of time steps dt and phase-space step 
dv that uniformly ensure a given precision of the numerical algorithm. 

Indeed, in both the time integration stage and the phase space integration 
stage, the integration method is an Euler scheme, either integrated in time or in 
the v variable. In these cases, it is well known that the precision of the method 
is bounded by the second derivative of the integrated variables multiplied by the 
integration step (the proof of this fact is close from the analysis performed in 
section [2]). More specifically, we have: 

(i). The numerical time integration consists of an Euler scheme approximating 
the solutions of the equations ([TJ on a compact subset of the phase space 
where \F(v) — w + 1\ < M. Moreover, the values of the variable (v,w) are 
contained in a compact set that can be chosen independent of the value of 9 



under assumption (A2) as already discussed. The second derivative of v(t) 
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and w(t) read: 

v "(t) = F'(v)(F(v) -w + I)-a{bv-w) + I'(t) 
w »(t) =a(b (F(v) - w + I))-a(bv- w) 

where x' denotes the derivative of x(t) with respect to time. We also note 
that the time-integration is performed on an even more reduced zone of 
the phase space. Indeed, the fact that \F(v) — w + I\ < M implies that 
F(v) < M + w u — I* defining an interval of values of v where the integration 
is performedm and on this compact set, we clearly have C\ > and Ci > 
such that \a{bv — w)\ < C\ and l-F'fv)! < C 2 . Therefore choosing dt < 
e/(max(C2M + C\ + ||/'||oo, \ab\M + C2)), we are ensured to have a precision 
of the order e on the evaluated value of v(t) and w(t). Moreover, besides the 
uniform bound on the error that allows choosing a time-step dt ensuring a 
precision of order e at any point of the phase plane where the time integration 
is performed, the analysis allows defining an optimal time step dt depending 
on the point (v, w) and ensuring the desired precision: 

dt(v, w) 



maxO"(f)|, \w"(t)\) 

This varying time step allows achieving the desired precision in a minimal 
number of operations. 

(ii). The numerical phase space integration consists of an Euler scheme approxi- 
mating the solutions of the equations |7f on a compact subset of the phase 
space. The second derivative of T(v) and W(v) with respect to v (denoted 
with a double dot) reads: 

■■ . , ab a(bv-W) aF'(v)(bv-W) 

Wlv) = 

yj F(v)-W + I (F(v)-W + I) 2 (F(v)-W + I) 2 

a(bv-W)(a(bv-W) + I'(T)) 



T(v) 



(F(v) -W + I) 3 
F'(v)(F(v) -W + I)-a(bv-W) + I'(T) 



(F(v) — W + I) 3 

and here again, we can find C3 > and C4 > such that: \a(bv — W)\ < 
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C 3 \F(v)-W + I\ and \F'(v)\ < C 4 \F(v) - W + I\ and if we denning dv such 
that: 

, ( r r ^ M+Cs + Cj \\r\\ooC 3 c A c 3 ^ ||/'||oo\ . 

we ensure that the error on W{y) and T{v) bounded by e. Here again, we 
are in position to choose an optimal integration step depending on the values 
of the variables (v, w) ensuring the desired precision: 

dv(v, w) = 



max(\W(v)\,\T(v)\) 



The evaluation of the values of C±, C 2 , C 3 , C 4 is quite difficult in the most 
general case. However, as mentioned, we can choose a varying integration steps dt 
and dv depending on the values of the variables that ensures exactly a precision 
bounded by e at each integration time, as done in the main text. 

It is interesting to note that the choice of M impacts the values of dt and dv 
and therefore the number of operations necessary to integrate the system for a 
fixed precision e. The larger M is chosen, the smaller dt is and the larger dv is, 
meaning that we emphasize the time integration (phase (i) of the algorithm) and 
give less importance to the phase space simulation (phase (ii) of the algorithm). 
Conversely, if M is small, the resulting dv will be smaller and dt larger, and the 
scheme will emphasize the phase-space integration. A balanced choice of M will 
allow an optimal simulation of the model by taking advantage of both integration 
methods, in the sense that it will perform the smaller number of operations for a 
given precision e. 
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